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ABSTRACT 


The  problem  of  estimating  the  arrival  times,  amplitudes,  and  phases  of  an 
unknown  number  of  signals  in  noise  is  treated.  The  signals  are  assumed  to 
have  a  common,  known  waveform.  A  Bayesian  model  using  a  Poisson  prior  for  the 
arrival  times  is  specified,  and  a  real  time  algorithm  for  computing  the 
posterior  mode  is  developed.  Alternatively,  the  procedure  may  be  looked  upon 
as  a  penalized  likelihood  estimator  with  a  penalty  term  which  is  a  generalized 
form  of  Akaike’s  Information  Criterion.  Simulation  results  are  presented 
which  show  that  this  approach  can  improve  over  classical,  linear  methods. 
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SIGNIFICANCE  AND  EXPLANATION 


The  problem  considered  here  is  the  extraction  of  an  unknown  number  of 
signals  from  noise  when  the  arrival  times#  amplitudes#  and  phases  are  unknown# 
but  that  the  signals  have  a  common#  known  waveform.  This  situation  arises  in 
sonar  and  radar  signal  processing.  The  signals  are  modelled  as  a  filtered, 
marked#  Poisson  process  and  it  is  proposed  that  the  maximum  posteriori 
estimate  be  used.  A  combinatorial  problem  arises  involving  the  sets  of 
possible  arrival  times#  which  has  plagued  many  similar  time  series  models.  An 
algorithm  is  proposed  for  computing  the  estimate  which  circumvents  thi3 
difficulty.  Simulation  results  show  the  procedure  can  improve  significantly 
on  classical  linear  procedures. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 
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ESTIMATION  OF  A  FILTERED  MARKED  POISSON  PROCESS  FROM  NOISY 
OBSERVATIONS  WITH  APPLICATIONS  TO  SIGNAL  PROCESSING 


Dennis  D.  Cox*1'2'3  and  John  E.  Ehrenberg**1 


1 »  Introduction 

Consider  the  problem  of  estimating  the  arrival  times,  amplitudes,  and 
phases  of  an  unknown  number  of  possibly  overlapping  signals  with  a  known  form 
of  finite  duration.  A  specific  example  of  such  a  problem  was  treated  by 
Ehrenberg,  Ewart,  and  Morris  (1978),  denoted  hereafter  as  EEM.  Briefly,  an 
acoustic  pulse  was  transmitted  underwater  to  a  receiver  some  distance  away. 

The  received  signal  actually  consisted  of  different  replicas  of  the 
transmitted  pulse,  each  arriving  at  a  different  time.  The  reason  for  this  is 
the  so  called  multipath  phenomenon  wherein  the  sound  is  transmitted  over 
several  paths  of  different  lengths.  The  observations  were  contaminated  by  the 
usual  ambient  noise  in  the  ocean  and  the  internal  noise  of  the  receiving 
apparatus . 

Similar  problems  arise  in  active  sonar  and  radar  wherein  one  sends  out  a 
pulse  of  energy  and  then  receives  the  echoes,  which  are  reflections  of  the 
pulse  off  of  objects  in  the  surrounding  environment.  The  echoes  have  a  known 
wave  form,  but  the  number  of  such  echoes,  their  arrival  times,  amplitudes,  and 
phases  must  be  estimated.  A  more  complete  description  of  this  problem  is 
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given  in  section  3.  A  mathematically  similar  problem  described  by  Barnard 
(1959)  concerned  a  production  process  which  undergoes  level  changes  of  unknown 
magnitudes  and  occurring  at  unknown  times.  The  problem  described  by  that 
author  was  different  from  the  one  considered  here  only  in  that  the  "signal"  (a 
level  change  or  step)  is  of  infinite  duration.  The  model  we  shall  propose 
below  is  a  simple  case  of  a  general  time  series  model  of  Abraham  and  Box 
(1979),  the  so  called  "aberrent  innovations"  model.  Those  authors  utilized  an 
approximation  wherein  the  number  of  signals  is  assumed  small  (<  2).  We  do 
not  allow  ourselves  this  luxury,  so  our  analysis  proceeds  somewhat 
differently. 

In  EEM,  a  maximum  likelihood  estimation  procedure  was  developed  under  the 
assumption  that  the  noise  is  white,  stationary,  and  Gaussian.  For  a  given  set 
of  arrival  times,  the  complex  amplitudes  of  the  associated  signals  were 
estimated  by  least  squares.  (The  magnitude  of  the  complex  amplitude 
corresponds  to  the  usual  amplitude,  and  the  argument  of  the  complex  amplitude 
corresponds  to  the  phase.)  For  a  given  number  of  signals,  all  possible 
combinations  of  arrival  times  were  examined,  and  the  one  with  highest 
likelihood  (minimal  residual  sum  of  squares)  was  chosen.  This  of  course 
involves  a  rather  heavy  computational  cost  for  two  or  more  signals  and  a  data 
record  of  even  moderate  length.  The  estimation  of  the  actual  number  of 
signals  presented  considerable  difficulty,  since  an  increase  in  the  number  of 
assumed  signals  results  in  an  increase  in  the  likelihood,  no  matter  what  the 
true  number  of  signals.  The  method  for  selecting  the  number  of  signals  used 
in  EEM  was  to  match  the  residual  sum  of  squares  to  the  noise  variance.  This 
method  seemed  to  work  quite  well. 

In  this  paper,  we  propose  a  Bayesian  modification  of  the  maximum 
likelihood  procedure  of  EEM,  and  an  efficient  algorithm  for  computing  the 


estimate.  In  Section  2,  a  stochastic  model  for  the  arrival  times  and  signals 


is  specified.  The  arrival  times  are  assumed  to  be  a  realization  of  a  Poisson 
point  process  (more  precisely,  a  discrete  time  approximation  thereto),  and  the 
amplitudes  of  the  signals  are  assumed  to  be  independent  complex  normal.  It  is 
then  proposed  that  we  take  as  estimate  that  collection  of  arrival  times  and 
associated  signal  parameters  which  maximize  the  posterior  likelihood.  The 
resulting  estimate  is  called  the  maximum  a  posteriori  (MAP)  estimate.  It  can 
also  be  described  as  a  penalized  maximum  likelihood  procedure.  The  likelihood 
of  EEM  is  penalized  in  two  ways.  Firstly,  there  is  a  penalty  on  the  signal 
amplitudes  from  the  independent  normal  assumption.  This  "shrinkage"  penalty 
appears  in  many  Bayesian  analyses.  A  second  penalty  results  from  the  Poisson 
process  prior  on  the  event  times.  It  is  this  penalty  which  determines  the 
estimated  number  of  events.  In  the  log  likelihood,  this  penalty  leads  to  a 
term  proportional  to  the  number  of  events  or  signals  (see  equation  (2.6) 
below) ,  and  is  thus  a  generalized  information  criterion  for  selecting  the 
number  of  events  (see  Akaike,  1974).  Such  a  penalty  for  the  likelihood  has 
been  proposed  for  a  different  problem  by  Bhansali  and  Downham  (1977). 

Bayesian  derivations  of  information  criterion  type  penalties  have  been  given 
in  other  contexts  by  Schwarz  (1978)  and  Akaike  (1978). 

A  main  contribution  of  this  article  is  the  description  of  a  relatively 
simple  algorithm  for  the  approximate  computation  of  the  MAP  estimate.  This  is 
described  in  Section  4.  The  simplicity  of  the  algorithm  relies  heavily  on  the 
assumption  that  the  signals  are  of  duration  which  is  short  compared  to  the 
overall  observation  record.  Hence,  the  change  in  likelihood  resulting  from 
inclusion  or  exclusion  of  a  particular  time  from  a  proposed  set  of  event  times 
can  be  calculated  using  only  the  observations  near  that  particular  time. 

Based  on  this  remark,  w<=  develop  a  hill  climbing  procedure  which  increases  the 


likelihood  with  each  iteration.  Furthermore,  because  the  estimate  depends  on 
the  observations  only  locally,  it  is  possible  to  do  the  computations  in  real 
time,  with  a  delay.  This  is  especially  important  for  signal  processing 
applications. 

One  of  the  main  difficulties  with  this  problem  comes  from  the  overlapping 
of  signals.  Indeed,  if  each  signal  lasts  only  one  discrete  time  unit,  then 
the  information  concerning  whether  a  particular  point  is  an  event  time  (and 
the  information  about  the  associated  signal  amplitude)  is  contained  in  the 
observation  itself  at  that  point.  However,  when  the  signals  are  of  a  duration 
of  two  or  more  time  units,  then  the  information  is  spread  out  over  the  whole 
observation  record,  because  it  is  necessary  to  account  for  signals  which  may 
overlap  with  a  possible  signal  at  the  time  of  interest,  and  to  account  for 
signals  which  may  overlap  with  those  which  may  overlap  with  the  signal  of 
interest,  etc.  In  terms  of  the  problem  of  EEM,  when  the  acoustic  pulses 
overlap,  the  resultant  may  have  greater  or  less  amplitude  than  the  components, 
depending  on  their  relative  arrival  times,  amplitudes,  and  phases.  This  is 
accounted  for  in  the  likelihood  equations,  but  determination  of  the  signal 
parameters  requires  solution  of  a  nondiagonal  linear  system  which  is  different 
for  each  possible  set  of  event  times.  Our  method  for  dealing  with  this 
difficulty  is  given  in  Section  4. 

In  Section  5,  we  present  performance  results  for  the  MAP  estimation 
procedure,  relying  primarily  on  simulations.  The  simulations  consisted  of  a 
signal  process  sample  path  generated  according  to  the  stochastic  model  for  the 
arrival  times,  including  complex  Gaussian  amplitudes,  with  Gaussian  white 
noise  added  in.  The  algorithm,  and  a  couple  of  competitors,  was  applied  to 
estimate  the  signal  process.  These  results  show  that  for  sufficiently  small 
values  of  the  event  intensity  parameter  (expected  number  of  signal  arrivals 
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per  unit  time),  the  MAP  estimate  will  significantly  outperform  the  classical 
linear  estimates.  For  larger  values  of  the  intensity,  there  will  be 
sufficient  overlapping  of  signals  that  the  underlying  process  will  lose  its 
impulsive  character  and  begin  to  look  like  an  ordinary  Gaussian  process  (see 
e.g.  Snyder,  1975)  for  which  a  linear  estimate  is  optimal.  The  most  important 
results  are  presented  in  Figures  1  and  2,  where  we  compare  the  predictive  root 
mean  squared  errors  of  our  procedure  and  the  best  linear  estimate.  Figure  3 
is  also  important  in  that  it  shows  that  the  error  is  relatively  insensitive  to 
misspecification  of  the  constant  in  the  generalized  information  jriterion  (by 
up  to  a  factor  of  8).  Other  results  presented  in  Section  5  indicate  how  well 
the  MAP  procedure  estimates  the  total  number  of  events  (Fig.  4),  the  number  or 
signals  whose  amplitude  exceeds  a  certain  threshold  (Figs.  5  and  6),  and  the 
arrival  time  of  the  largest  amplitude  signal  (Fig.7).  it  is  only  with  respect 
to  this  latter  criterion  that  the  classical  methods  appear  to  offer  any 
serious  competition  for  the  MAP  procedure. 
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2.  Model  Description  and  Likelihoods 


In  this  section  we  specify  a  discrete  time  stochastic  model  for  the 

signal  process  described  in  the  introduction.  This  model  will  be  related  to  a 

continuous  time  model  for  the  active  sonar  problem  in  the  next  section. 

Because  this  application  involves  sinusoid  1  signals,  it  will  be  expedient  to 

use  complex  valued  stochastic  processes.  An  asterisk  *  will  denote  the 

operation  (on  vectors  or  matrices)  of  taking  the  complex  conjugate 

2 

transpose.  Given  a  complex  y  and  a  >  0,  by  a  complex  Gaussian  random 

2  2 

variable  Z  with  mean  y  and  variance  o  (written  Z  *CK(p,c  )),  we  mean 

that  the  real  and  imaginary  parts  of  Z  (written  He  Z  and  Im  Z)  are 

2 

independent  real  valued  random  variables  with  distributions  N(Re  y,o  )  and 
2 

N (Im  y,cr  ),  respectively. 

The  observed  vector  is  a  (complex)  N  vector  which  is  assumed  to 
satisfy  the  regression  equation 

(2.1)  £  =  Ka  +  £' 

where  the  noise  vector  e  has  independent  and  identically  distributed 

2 

(abbreviated  i.i.d. )  components  with  the  CN(0,o  )  distribution,  W  denotes 
an  N  X  N  Toeplitz  matrix  with 

(2.2)  Wjk  =  w(j-k) 

for  some  complex  function  w  defined  on  the  integers  (referred  to  as  a 
waveform) ,  and  a  is  an  N-vector  of  complex  amplitudes  whose  components  are 
i.i.d.  with  the  mixture  distribution 

(2.3)  ak  ~  p  CN(0 , p2 )  +  (1-p)  SQ  . 
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Here,  0  <  p  <  1,  p  >  0,  and,  6^  denotes  a  unit  point  mass  at  0.  Thus, 

2 

with  probability  p,  a^  is  chosen  according  to  the  CN(o,p  )  distribution, 
while  with  probability  t  -  p,  ak  is  0.  We  further  assume  that  for  some 
positive  integer  m,  the  waveform  w  satisfies 


(2.4) 


w(k)  *  0  ii  k  <  0  or  k  >  m. 


so  that  W  is  an  (m+1)  -  banded  matrix.  We  assume  that  the  parameters 
2  2 

a  ,  p,  p  ,  and  w( • )  are  known,  and  consider  the  problem  of  estimating  a 
given  the  observation 

Under  the  assumptions  set  forth  above,  we  have  a  Bayesian  regression 
problem,  and  the  posterior  distribution  of  a  given  £  can  be  obtained  in 
the  usual  way.  In  order  to  give  a  simple  description,  define  the  (random)  set 


of  integers 


K  “  {k:  1  <  k  <  N,  a^  *  0}. 


In  other  words,  K  is  the  set  of  event  times,  i.e.  the  times  k  when  ak  is 

2 

chosen  from  the  CN(0,p  )  distribution,  at  least  with  probability  one.  We 


also  put 


|K  |  =  number  of  elements  (cardinality)  of  K 


and  write  K  in  the  form 


and  put 


K  =  {k( 1 ) ,  k(2) , . .. ,k( |K |) } 


k(1)  <  k(2)  <  ...  <  k(|K|), 


a  (a(k( 1 ) ) ,  ...,  a(k(  |K  | ) ) )' . 


We  shall  frequently  pass  from  subscript  notation  to  functional  notation  as  in 
this  last  equation  wherein  we  wrote  a(k)  for  ak.  We  also  define  the 
H  x  |K|  matrix  W^,  which  is  the  submatrix  of  W  obtained  by  deleting  the 
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columns  whose  indices  are  not  in  K  .  Hence,  the  equation 


Wd  *  W[/  ft  1/ 

holds • 

The  posterior  density  for  K  ,  a has  the  form 


(2.5) 


f(K,aK|£)  =  c(£)  •  Cp/d-p))^1  •  {2np2)~|K|  • 


-2  *  *  i  *  -2  -2  * 

exp[a  Re  a^,  Wj,  (p  %|+  0  vk  * 


IK  I 


K  K'  n 


2  2 

Here,  c(j£)  depends  only  on  jr  (and  o  ,  p,  p  ,  and  w)  but  not  on  K  and 
a^,  and  1^  is  an  |n|  x  |n|  identity. 

Given  the  posterior  density  (2.5),  it  is  possible  at  least  in  principle 
to  obtain  an  optimal  Bayes  estimate  of  the  amplitude  vector  a,  such  as  the 
posterior  mean.  This  will  involve  a  complicated  integration.  For  each 
K  C  {1,2,...,N},  it  is  possible  (in  principle)  to  integrate  out  a^,.  This 
involves  computing  the  inverse  and  determinant  of  a  matrix  which  depends 
on  K  .  However,  there  are  2N  values  of  K,  so  that  an  exact  procedure  is 
impractical  for  even  moderate  sample  sizes  N.  Under  certain  circumstances, 
the  problem  is  tractable  (e.g.  when  signal  length  is  1,  i.e.  m  =  0),  and  there 
are  various  approximations  that  would  work  sometimes.  For  example,  if  p  is 
sufficiently  small,  then  only  values  of  K  for  which  |K|  <1  would 
contribute  to  the  integral.  However,  we  are  interested  in  the  general  case 
where  such  approximate  methods  would  not  be  accurate. 

One  estimate  which  is  usually  easier  to  compute  than  the  posterior  mean 
is  the  posterior  mode,  also  known  as  the  maximum  a  posteriori  (MAP) 

A 

estimate.  In  this  situation,  the  MAP  estimate  K ,  are  those  values  of 

TV 


K, which  maximize  the  posterior  likelihood  (2*5),  or  equivalently,  minimize 

(2.61  L(K,aKlX>  ■  4  4-  +  1/2 %’  h  +  ll:|a 

=  -Re  a(K)*}£  +  V2a(K)*Q  a(K)  +  IK  I  a. 

Here  we  have  introduced  the  not\tions 

jj  =  w*  £  e  cN 

l  A' I 

=  w*  x  e  c'M 

(2.7)  Q  *  (a/p)2iN  +  w*w 

2r(0/pl\l+WKWK 

a  =  -a2  log  {p/(2irp2  [1-p] )  }. 

We  have  also  written  a(K)  for  a  e  CN  to  make  explicit  the  set  K  of 
indices  corresponding  to  nonzero  entries. 

Now,  for  a  given  set  K  £  {1/2/.../N},  we  may  obtain  the  MAP  estimate  of 
a,  denoted  4(K)  e  (or  €  C1*1  if  the  zero  components  are  deleted)  by 
solving  a  |K  |  x  |K  I  linear  system.  The  result  is 

(2*8)  k  =  gk  k 

and  then  the  negative  log  likelihood  becomes 

uk,  &K  |  X)  =  - V2a *  H/c+  Ik i o. 

We  are  then  obliged  to  search  over  all  2W  values  of  K  to  find  a  K 
which  minimizes  the  last  expression,  leading  us  back  to  the  same  dilemma  as 
presented  by  the  posterior  mean  computation.  Since  p  *  0,  it  is  reasonable 

A 

to  suppose  that  o  >  0,  whereas  if  o  <  0,  we  simply  have  K  =  (1,2,...,N}. 
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In  Section  4  below,  we  present  a  hill  clitrbing  algorithm  which  attempts  to 

A 

find  K,  or  an  approximation  thereto. 

We  would  like  to  now  briefly  point  out  some  of  the  salient  features  of 

the  MAP  estimation  procedure.  Firstly,  note  that  the  problem  of  determining 

for  a  given  K  (the  solution  given  in  (2.8))  is  a  standard  regression 

problem,  although  our  use  of  a  Gaussian  prior  introduces  some  shrinkage  (the 
2 

(0/p)  Ij/  terms  in  (2.7)).  A  non-Bayesian  may  be  content  to  consider  this  as 
a  ridge  regression  type  estimate.  There  is  now  a  considerable  body  of  litera¬ 
ture  on  how  to  estimate  the  ridge  parameter  (or  the  hyperparameters  in  the 

2 

Bayesian  model),  and  such  estimates  could  be  transformed  into  estimates  of  p  . 

Within  this  regression  context,  the  estimation  of  the  set  K  of  event 
times  translates  into  a  model  selection  problem,  wherein  K  corresponds  to 
the  set  of  variables  in  the  model.  One  then  recognizes  (2.6)  as  a  penalized 
(posterior)  log  likelihood,  where  the  penalty  term  lK|a  corresponds  to  a 
generalized  version  of  AkaJke's  Information  Criterion  (AIC),  as  described  in 
the  introduction 

We  will  now  describe  two  other  procedures  which  could  be  applied  to  this 
problem,  which  we  refer  to  as  the  matched  filter  estimate  (abbreviated  MFE) 
and  best  linear  estimate  (BLE).  Still  other  methods  are  described  in  EEM.  If 
one  believed  that  there  was  an  event  at  time  k  whose  effect  did  not  overlap 
with  any  others,  then  an  unbiased  least  squares  estimate  of  the  amplitude  is 

\ ;  V[5kx  -  l0/p)2)- 

In  one  form  or  another,  such  estimates  are  widely  used  in  signal  processing, 

and  will  be  referred  to  as  matched  filter  estimates.  Our  MAP  procedure  is  a 

2 

modification  of  this  matched  filter  in  that  the  shrinkage  term  ( a/p)  is 
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included  in  the  denominator,  and  an  effort  is  made  to  take  into  account  the 
overlaps  of  the  various  effects  (which  leads  to  the  system  of  equations  whose 
solution  is  given  in  (2.8)).  Indeed,  if  there  is  no  overlap  (i.e.  if  m  =  0). 
then  the  MAP  estimate  is  essentially  the  MF  estimate  for  those  times  where  the 
MFE  is  above  a  suitable  threshold  determined  by  the  information  penalty. 

Under  these  circumstances,  the  MAP  procedure  just  amounts  to  simultaneous 
signal  detection  and  estimation.  However,  when  overlapping  occurs,  the  MFE 
becomes  quite  confused  and  will  indicate  the  presence  of  large  amplitude 
signals  when  sufficiently  many  small  signals  overlap  and  add  more  or  less  in 
phase.  The  MAP  procedure  provides  a  means  for  correcting  this  defect. 

Another  approach  would  be  to  use  the  linear  estimate  which  is  optimal  for 
the  Bayesian  model  (over  all  linear  estimates)  in  the  sense  of  minimizing  mean 
squared  error.  This  estimate  is  the  minimizer  over  a  6  C?1  of  the  quadratic 
function 

(2.9)  1/(2pp2)  a*a  +  1/(2o2)(jr  -  Wa)*(£  -  Wa ) . 

We  will  refer  to  this  estimate  as  the  BLE. 

Both  the  MF  and  BL  estimates  have  the  difficulty  that  they  don't  take 
into  account  the  impulsive  nature  of  the  stochastic  model  for  a.  Both 
procedures  will  almost  never  estimate  aj  =  0  for  some  j,  and  hence  do  not 
immediately  give  estimates  for  the  event  times.  Also,  both  methods  have 
difficulty  separating  a  large  number  of  small  overlapping  signals  from  a 
single  large  signal,  and  they  will  give  large  amplitude  estimates  at  times 
near  a  large  signal  rather  than  just  at  its  single  event  time.  The  MAP 
estimate,  however,  does  not  suffer  from  these  problems. 
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In  this  section  we  briefly  describe  a  class  of  real  problems  which 
motivated  th  ,>  current  study,  and  show  how  to  link  them  up  with  the  model  of 
the  previous  section.  The  use  of  sonar  and  radar  models  involving  filtered 
Poisson  processes  is  not  new.  Middleton  (1979)  has  developed  a  very  general 
model  incorporating  such  a  noise  component  and  has  written  at  length  on  the 
various  estimation  problems  for  such  models.  See  also  Ehrenberg  (1976).  In 
active  sonar,  a  narrowband  sonic  pulse  of  known  frequency  fQ  sec.  and 
with  a  known  envelope  function  wQ(t)  is  transmitted,  and  the  echoes  formed 
by  reflections  of  the  sonic  pulse  off  of  solid  objects  are  received.  The 
signal  process  which  results  from  many  such  echoes  is  commonly  called 
reverberation  in  the  sonar  literature  and  clutter  in  the  radar  literature. 

The  echoes  have  unknown  return  times  (depending  on  the  distance  of  the 
reflecting  object)  and  amplitudes  (depending  on  the  size  and  distance  of  the 
object),  but,  if  the  reflecting  objects  are  stationary  so  that  there  is  no 
Doppler  shift,  the  echoes  are  of  the  same  frequency  and  shape  as  the 
transmitted  pulse.  Therefore,  the  received  echoes  as  a  function  of  time 
constitute  a  "signal"  of  the  form 

S(t)  =  Re  l  A(j)  wQ(t  -  T(  j ) )  exp[2irifQ  (t  -  t(  j ) )  ] 

where  A(j)  and  x(j)  are  the  complex  amplitude  and  arrival  time  of  the  j^ 
echo.  In  practice,  there  is  always  noise  contaminating  the  signal,  so  that 


the  observation  model  is 


where  N(*)  denotes  the  noise,  which  we  assume  is  white  Gaussian  noise  with 
two  sided  power  spectral  density  y  >  0. 

The  initial  stage  of  the  signal  processing  involves  reduction  of  the 
continuous  time  data  to  the  discrete  time  quadrature  demodulator  statistics 
given  by 


k  T/N 

yv  “  m  /  y  (t)  exp  [-2itif  t]  dt,  for  1  <  k  <  N. 

(k-1 )T/U 

It  is  assumed  that  fgT/N  is  an  integer,  and  also  that  wQ  is  nearly 
constant  over  a  time  interval  of  length  T/N.  Under  the  white  noise 
assumption,  the  noise  vector  e  is  given  by 

kT/N 

£;  =  -  /  N(t)  exp  [~2itif  t]  dt 

(k-1 )T/N  0 

2 

where  these  components  are  i.i.d.  CN(0,o  )  distribution,  with 

(3.1 )  a2  =  YN/T. 

Now  suppose  that  the  echo  arrival  times  x( 1 ) ,  t(2),...  are  the  event 
times  of  a  homogeneous  Poisson  process  with  constant  intensity  X,  and  thut 

(3.2)  p  =  X  T/N 

is  sufficiently  small  that  we  may  assume  that  at  most  one  echo  arrives  during 
one  of  the  Fourier  processing  intervals  of  duration  T/N.  Further  assume  that 
Wq{»)  is  sufficiently  smooth  so  as  to  be  nearly  constant  over  each  Fourier 
processing  interval  and  then  set 
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(3.3) 


w(k)  =  wQ(k  N/T) 


Finally,  suppose  that  fQ  T/N  >>  1  so  that  the  phase  of  the  arriving  signal 

is  essentially  random,  and  that  the  magnitude  |A(j)|  of  A(j)  has  the 

2 

Rayleigh  distribution  with  parameter  p  (see  Middleton,  1960  for 

definitions),  and  hence  the  complex  amplitude  of  the  arriving  signal  will  have 
2 

the  CN(0,p  )  distribution.  This  discussion,  together  with  equations  (3.1) 

through  (3.3),  specifies  how  to  fit  this  sonar  problem  into  framework  of  the 

model  defined  in  Section  2.  We  note  that  the  problem  treated  in  EEM  may  also 

be  fit  in  this  framework  in  an  obvious  way,  if  the  amplitudes  are  Rayleigh. 

A  possibly  better  model  of  the  active  f  Dnar  problem  would  be  obtained  by 

assuming  that  the  solid  scatterers  which  produce  the  echoes  have  a  homogeneous 

spatial  Poisson  distribution  with  parameter  a  per  unit  of  volume  and  that  at 

a  unit  distance  from  the  transmitter,  their  echo  amplitude  would  be 
2 

CN(0,p^  ).  Also  suppose,  for  the  sake  of  simplicity,  that  the  sonar  transmit 
beam  is  conical,  subtending  a  solid  angle  of  measure  <o.  Then  the  echo 
arrival  time  process  is  a  nonhomogeneous  Poisson  process  with  rate  function 

3  2 

X(t)  =  a  toe  t  , 

c  being  sound  velocity,  and  the  complex  amplitude  of  an  echo  arriving  at  time 
2 

x  is  CN(0,p  (t),  where 


P(t)  =  P1 /( c t)  . 


These  complications  can  be  easily  incorporated  in  the  model  of  section  2,  and 
will  change  the  log  likelihood  (2.6).  For  example,  the  term  |fC  |  a  will  be 
replaced  by 
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where 


l  atk) 
k  e  K 


a(k)  =*  -log  [{X(kT/N)T/N}/{2irp  (kT/N)(1-X(kT/N)T/N)}]. 


We  maintain  that  the  necessary  modifications  are  more  notational  than 
substantial,  and  also  that  an  understanding  of  the  MAP  estimation  problem 
posed  in  Section  2  will  immediately  lead  to  an  understanding  of  the 
corresponding  MAP  estimation  procedure  for  this  more  complicated  problem. 


4«  Recursive  Computation  of  the  MAP  Estimate 

In  this  section,  we  describe  an  algorithm  for  approximate  computation  of 
the  HAP  estimate  described  in  Section  2.  This  algorithm  starts  with  an 
initial  estimate  of  K  and  the  corresponding  given  by  (2.8)  (e.g. 

initially  K  =  fS,  in  which  case  no  5^,  is  required),  and  produces  a  final 
estimate  of  K,  together  with  the  corresponding  which  has  higher 

posterior  likelihood.  One  of  the  main  features  of  the  procedure  is  that  it 
can  be  implemented  in  such  a  way  that  only  one  pass  through  the  data  is 
required,  and  such  that  at  any  time  only  a  small  span  of  the  data  is  being 
processed.  This  is  significant  to  statisticians  with  large  data  sets,  and  to 
engineers  who  are  interested  in  real  time  processing. 

We  now  give  an  overview  of  the  algorithm  via  an  hierarchal  description, 
proceeding  from  macro  to  micro.  The  interested  reader  may  consult  the  Fortran 
program  segment  listed  in  Appendix  A. 

(i)  There  are  several  (say  n  )  applications  of  the  posterior  likelihood 
increasing  procedure  described  in  stage  (ii).  These  can  be  working 
simultaneously  on  adjacent  sections  of  rhe  data  stream. 

(ii)  Local  perturbations  in  the  set  K  of  estimated  event  times  are 

successively  examined.  With  each  such  perturbation,  the  updated  amplitude 
vector  estimate  and  the  log  likelihood  increment  are  computed  as 

described  in  (iii)  and  (iv).  The  simplest  example  of  such  a  perturbation  is, 
for  a  particular  time  epoch  j ,  to  augment  the  current  K  with  j  if 

j  &  K  ,  or  to  delete  j  from  K  if  j  e  K,  Afterwards,  the  negative  log 
likelihood  is  checked  for  a  decrease.  If  there  is  one,  the  K  is  updated, 
and  otherwise  it  is  left  unchanged.  The  processor  would  then  proceed  to  the 
(j  +  1)st  epoch  to  repeat  the  procedure  there,  and  so  forth. 
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(iii)  Each  time  K  is  perturbed,  say  from  KQ  to  f(1 ,  it  is  necessary 
to  update  &£,  from  a(KQ)  to  &(£.,),  which  involves  solving  a  new  linear 
system  (see  (2.8)).  A  method  for  doing  this  is  described  in  (iv).  The 
increment  in  the  negative  log  likelihood  is  evaluated  by 

AL^-V^S^)  -  4(K0))*  E+  (1K1  i  -  |K0I>  a. 

If  AL  <  0,  then  Kq  is  replaced  by  and  otherwise  Kq  is  left 

unmodified. 

(iv)  Given  KQ,  £(KQ),  and  K^,  it  is  required  to  compute  5(K  1 ) . 

Since  K i  is  obtained  by  a  local  perturbation  on  Kq,  4(fC.  )  should  be 
obtainable  by  locally  updating  S(| <Q).  This  can  be  accomplished  by  a  few 
applications  of  a  Gauss-Seidel  recursion.  A  single  iteration  involves 
replacing  the  kth  component  of  £  ^  by  the  updated  value 

Yk  =  (UfC,k  "  Qf(,kj  \',j)/Qf<,kk 

Since  is  2m  +  1  banded,  there  are  really  only  2:n  terms  in  the 

summation. 

We  now  indicate  the  span  of  the  data  stream  which  is  involved  in  the 
updating  at  any  one  time.  Each  Gauss-Seidel  iteration  involves  (2m  +  1) 
adjacent  data  points,  and  if  M  such  iterations  are  performed,  then  2mM  +  1 
adjacent  data  points  are  included.  Note  that  if  we  are  updating  £  as  a 
result  of  a  perturbation  of  K  at  j,  this  entire  updating  procedure  must  be 
completed  before  moving  on  to  the  perturbation  at  j  +  1.  However,  the  Gauss- 
Seidel  iterations  can  be  condensed  into  a  single  matrix  multiplication  of 
the  2mM  +  1  data  points  by  a  square  matrix.  Now  if  n  of  the  overall 


likelihood  increasing  steps  in  (i)  are  applied,  then  a  span  of  n(2mM  +  1) 
data  points  are  actually  being  processed  at  one  time.  Our  simulation 
experience  (results  of  which  are  reported  in  the  next  section)  indicate 
n  *  2  and  M  »  3  are  usually  sufficient. 

We  now  indicate  two  modifications  of  the  basic  algorithm  which  were  used 
for  the  simulations  reported  in  the  next  section.  One  can  expand  the  local 
perturbation  in  the  set  K  as  follows.  For  a  given  time  epoch  j  of 
interest,  one  can  examine  all  8  combinations  of  inclusion/exclusion  of  the 
epochs  j ,  j  +  1 ,  and  j  +  2  from  the  set  K  ,  compute  the  likelihood  change 
for  each,  and  modify  K  accordingly.  Then  at  the  next  step,  look  at  j  +  1 , 
j  +  2,  and  j  +  3,  etc.  This  procedure  can  be  extended  to  looking  at  blocks 
of  times  of  length  4,  5,  etc.  Also,  it  is  not  necessary  to  use  successive 
blocks  of  epochs,  one  could  say  skip  from  j  to  j  +  2,  as  long  as  there  was 
substantial  overlap  of  the  blocks.  However,  we  stayed  with  the  simple  version 
of  successive  blocks  of  3.  Another  modification  is  to  run  through  the  data 
with  the  Gauss  Seidel  smoothing  procedure,  but  keeping  the  current  K 
fixed.  This  serves  to  give  better  accuracy  in  solving  the  linear  system  of 
equations  corresponding  to  K  . 
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5 .  Simulation  Results 

In  this  section,  we  present  simulation  results  for  the  MAP  estimate.  The 
purpose  of  these  simulations  is  to  show  that  the  MAP  procedure  can  be 
implemented  with  a  small  number  of  each  cf  the  basic  steps  and  still  perform 
quite  well  in  comparison  with  other  procedures.  The  good  performance  should 
come  as  no  suprise,  given  the  results  presented  in  EEM.  What  is  significant 
is  that  our  algorithm  achieves  this  level  of  performance.  The  MAP  estimate 
was  compared  wherever  possible  with  the  MF  and  BL  estimates,  and  only  for  one 
performance  measure  (estimating  the  maximum  amplitude  event  time)  were  they 
serious  competitors.  There  are  four  different  performance  criteria  which  we 
consider,  which  are  estimation  of  the  total  signal  (the  filtered  point 
process,  or  the  superposition  of  all  the  individual  signals),  the  total  number 
of.  signals,  the  number  of  signals  with  amplitudes  larger  than  a  certain 
threshold,  and  the  event  time  of  the  maximum  amplitude  signal.  We  also 
considered  the  performance  of  the  MAP  estimate  when  the  information  penalty 
factor  a  (see  (2.7))  is  misspecif ied,  and  found  that  it  was  surprisingly 
insensitive  to  this  misspecif ication.  This  is  also  apparent  in  the  simulation 
results  of  Bhansali  and  Downham  (1977).  Perhaps  this  is  one  of  the  reasons 
AIC  has  worked  so  well  in  practice. 

We  now  specify  some  of  the  parameters  of  the  Monte  Carlo  experiment.  All 
points  plotted  in  the  accompanying  figures  represent  an  average  of  25  trials 
of  N  =  100  observations  each.  In  all  cases,  the  waveform  w  was  unity,  i.e. 

(5.1)  w(j)  =  1  for  0  <  j  <m. 

The  signal  length  (m  +  1 )  was  either  4,  7,  or  10.  We  use  L  to  denote 

2 

(m  +  1)  for  convenience.  The  noise  variance  o  was  fixed  at  1.  The  event 
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probability  p  was  always  chosen  as  a  power  of  V? ,  either  .25,  .125, 

.0625,  or  .03125,  corresponding  to  values  of  2,  3,  4,  or  5  for 

-log2  p.  Note  that  the  expected  number  of  signals  being  received  at  any 

2 

particular  time  is  just  Lp.  The  amplitude  variance  p  was  chosen  so  that 
the  quantity 

2 

c  =  pp 

took  on  one  of  the  values  .25,  .5,  1,  or  2.  Hence,  as  p  decreases 

(-log2  p  increases),  the  arriving  signals  become  fewer  and  larger  in 

amplitude.  Note  that  the  process  covariance  is  determined  by  c,  L,  and 
2 

0  ,  so  that  specifying  c  was  convenient  for  the  BLE.  This  also  means  that 
the  signal  to  noise  ratio  is  the  same  for  the  simulations  at  the  same  value 
of  c.  In  all  cases,  n  =  3  passes  with  the  basic  perturbation  algorithm 
were  made  (step  (i)  of  the  previous  section),  and  M  -  3  Gauss  Seidel 
iterations  (step  ( iv) )  were  performed  at  each  step.  Also,  there  was  one  final 
pass  through  the  data  with  a  Gauss-Seidel  iteration,  as  described  in  the  last 
section. 

The  first  performance  criterion  we  examine  is  the  integrated  mean  squared 
error  (IMSE)  in  the  estimate  of  total  signal  process,  where  the  total  signal 
process  is  taken  to  be  the  sum  of  all  the  individual  signals  (received 
echoes).  More  formally 

IMSE  =  E  ||  W(S  -  a) | |2, 

where  ||*||  denotes  the  usual  norm  on  C*1.  T>._i  ordinates  in  Figures  1 

through  3  are  the  square  root  of  the  IMSE  (referred  to  as  just  the  RMS  error), 
divided  by  the  total  RMS  signal  power,  which  is  given  by 
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2  1/, 

(RMS  signal  power)  *  tNpL(2p  )]  2=>  [2NcL]  2  . 


2 

Note  that  the  totrl  noise  power  is  just  2No  .  The  estimates  of  the 
normalized  RMS  error,  each  based  on  a  sample  of  25  simulation  trials,  are 
reasonably  accurate.  In  no  case  did  the  (estimated)  standard  error  exceed 
.025,  and  was  usually  less  than  .015.  Three  points  in  Figure  1  were 
replicated,  and  one  notes  how  close  the  values  were  in  each  case. 

On  Figures  1  and  2,  we  have  also  shown  the  normalized  steady  state  RMS 
error  for  the  BLE.  Using  standard  spectral  methods  as  in  Koopmans  (1974),  the 
steady  state  MSE  per  observation  is  given  by  (see  Appendix  B) 


r  /  — -1  a> 

*  -x  c  V  +  |w(X)|2 


where 


1  w( X) | 2  =  sin2  (LX/2)/sin2  ( X/2) 


This  will  in  fact  be  a  lower  bound  for  the  exact  IMSE  for  the  BLE  for  a  finite 

data  record  because  of  boundary  effects,  and  the  accuracy  will  decrease  with 

increasing  L  or  c.  However,  we  believe  these  numbers  are  reliable.  The 

main  conclusion  to  be  drawn  from  Figures  1  and  2  is  that  for  p  <  ,08,  the 

MAP  estimate  does  better  than  the  BLE,  in  terms  of  IMSE.  For  L  =  10  and 

p  =  .08,  one  expects  .8  signals  being  received  at  any  one  time  (total  of  8 

expected  events  in  the  whole  record) ,  so  there  is  considerable  overlapping  of 

signals.  Also,  on  these  two  figures  are  shown  some  limiting  IMSE  values  for 

2 

the  MAP  estimate  as  p  +  0  (holding  pp  =  c  constant).  Under  these 
conditions,  there  will  be  essentially  no  overlapping  effects  and  the  event 
time  set  K  will  be  perfectly  estimated.  It  can  then  be  shown  (Appendix  C) 
that 
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f 


IKSE(p+0)  -  2NpL02/(02/p2  +  L). 

In  any  event/  this  is  always  a  lower  bound  for  the  MAP  IMSE,  with  decreasing 
accuracy  for  increasing  L. 

Figure  3  shows  a  plot  for  the  normalized  RMS  error  when  the  penalty  a 
is  mismatched.  The  "correct"  value/  denoted  just  a,  is  given  in  (2.7),  but 
this  is  not  necessarily  optimal.  Denoting  by  a  the  value  which  was  in  fact 
used/  we  set  a  ■  a  in  all  simulations  except  for  those  reported  in  Figures  3 
and  6.  It  is  clear  from  Figure  3  that  a  is  very  nearly  the  optimal  value 
(for  minimizing  IMSE)  and  furthermore  that  the  IMSE  of  the  estimate  is  quite 
insensitive  to  the  actual  a  used. 

Our  second  performance  criterion  is  the  error  in  estimating  the  total 
number  of  individual  signals.  In  Figure  4,  the  ordinates  are  the  RMS  error  in 
the  estimate  of  the  total  number  of  signals  in  the  simulated  data  record  of 
100  observations,  normalized  by  dividing  by  (Np{1  -  p))^2,  the  standard 
deviation  of  the  number  of  signals.  Note  that  this  latter  quantity  would  be 
the  RMS  error  if  one  used  the  constant  Np  as  the  estimate  of  the  number  of 
signals.  There  is  considerable  sampling  variability  in  the  results  of  Figure 
4,  as  is  clearly  seen  by  comparing  the  replicated  points.  However,  for 
p  <  . 1 ,  it  appears  that  the  MAP  estimate  does  do  some  good.  In  simulations 
with  mismatched  a,  the  performance  with  respect  to  this  indicator  degraded 
quite  severely.  The  reason  for  this  is  of  course  the  fact  that  there  are  many 
signals  of  small  magnitude  which  are  difficult  to  accurately  detect,  and  if 
a  is  set  very  low  to  detect  them,  then  many  false  detections  are  made.  Also 
note  that  there  is  no  simple  way  to  use  the  MFE  or  BLE  to  estimate  the  total 
number  of  signals. 
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Figures  5  and  6  present  the  RMS  error  in  the  MAP  estimate  of  the  number 

of  signals  for  which  the  amplitude  i«jj  exceeded  p.  The  estimate  was 

obtained  by  counting  the  number  of  MAP  amplitude  estimates  exceeding  p  in 

magnitude.  For  a  given  signal,  the  probability  its  amplitude  will  exceed  p 
2 

is  just  P(x2  >  1]  *  e-ty2  “  .6065  ,  so  . 6065Np  such  signals  are  expected, 
and  the  standard  deviation  of  the  number  of  such  signals  is 

[.6065Np  (i  -  .eoesp)]1^ . 

The  KMS  errors  o.!  the  MAP  estimate  were  divided  by  this  last  quantity  for 
normalization.  In  Figure  5,  we  see  that  the  MAP  estimate  is  not  significantly 
better  than  the  constant  estimate  of  .6065Np  for  p  =*  .25,  but  exhibits 
sharp  improvement  as  p  decreases.  Compare  the  slope  in  Figure  5  with  that 
in  Figure  4.  We  also  note  from  Figure  6  that,  at  least  for  p  =  .0(525,  the 
performance  is  relatively  stable  under  misspecification  of  the  information 
penalty  fact  a  by  a  factor  of  2.  In  that  figure,  the  point  p  =  .125, 
a/a  =  .5  was  replicated  because  the  first  value  (with  error  /st.  dev.  value 
of  .9)  indicated  a  discontinuity  in  going  from  "a/a  =  .5  to  1.  However,  the 
second  value  did  not  clear  up  matters  very  much.  These  results  of  Figure  6 
are  important  to  those  who  are  interested  in  just  estimating  the  number  of 
signals,  or  the  event  probability  p,  since  p  only  enters  in  the  procedure 
through  a.  It  is  of  course  possible  to  use  the  BL  and  MF  estimates  to 
estimate  the  number  of  signals  with  amplitudes  exceeding  p,  and  we  did  in 
fact  simulate  these.  The  BLE  was  much  better  than  the  MFE  in  this  regard,  but 
was  always  much  worse  than  the  MAP  estimate  with  RMS  errors  two  or  more  times 
those  of  the  MAP  estimate. 


-23- 


Our  final  performance  criterion  concerns  the  estimation  of  the  arrival 
time  of  the  largest  amplitude  signal,  abbreviated  MAET  for  maximum  amplitude 
event  time.  For  each  type  of  estimation  procedure,  the  MAET  was  estimated  by 
the  time  of  the  largest  estimated  amplitude.  Some  results  are  presented  in 
Figure  7,  where  the  ordinate  represents  the  fraction  of  25  simulation  runs 
that  the  procedure  picked  the  correct  MAET.  We  should  note  that  with  very  few 
exceptions,  each  of  the  procedures  chose  exactly  the  correct  time,  or  was  off 
by  three  or  more  time  periods.  Our  expectation  was  that  the  MFE  and  BLE  would 
be  frequently  "fooled"  by  overlapping  signals  which  add  to  give  a  large 
amplitude,  whereas  the  MAP  would  separate  out  the  individual  signals. 

However,  the  results  seem  to  indicate  rough  equality  of  the  three  different 
estimates.  We  do  note  that  there  is  an  apparent  trend  for  the  MAP  estimate  to 
perform  better  than  the  other  two  for  small  values  of  p.  Unfortunately,  it 
would  take  a  large  number  of  simulation  trials  to  validate  this. 
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6«  Summary 

The  statistical  problem  treated  here  has  been  investigated  by  other 
authors,  but  there  has  been  relatively  little  progress  for  moderate  values  of 
the  event  time  intensity  parameter  p.  This  is  because  of  the  apparent  need 
to  examine  at  least  a  large  fraction  of-  the  collection  of  all  sets  of  event 
times.  We  have  proposed  an  estimation  procedure  which  avoids  this  difficulty, 
and  in  fact  can  be  implemented  in  a  real  time  algorithm.  Simulation  results 
indicate  that  the  procedure  can  do  much  better  in  many  respects  than  classical 
methods.  We  hope  that  this  work  will  stimulate  others  to  improve  upon  the 
basic  algorithm,  and  to  develop  a  better  understanding  of  its  statistical 

properties. 
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Appendix  At  An  ASCII  Fortran  program  segment  for  implementing  the  MAP 
estimation  algorithm.  The  complex  array  A  stores  the  current  set  of 
amplitude  estimates,  with  many  elements  equal  to  0,  and  the  array  B  is 
used  for  scratch  work  to  determine  the  effect  of  a  perturbation  in  K .  The 
internal  subroutine  MAPXNC  determines  the  effect  of  the  change  in  K,  using 
the  Gauss-Seidel  computations  in  CINC.  The  vector  p  is  in  RMU,  and  the 
matrix  Q  is  stored  in  a  1  dimensional  array  since  it  is  Toeplitz  and 
Hermitian.  The  outer  loop  (DO  20)  is  to  implement  (i),  while  the  big 
IF. . .THEN. . .ELSE  implements  (ii).  The  change  in  negative  log  posterior 
likelihood  (DLIK)  is  computed  by  LIKCHG.  If  the  change  is  negative,  the 
A  is  replaced  by  B  with  BTOA,  and  otherwise  B  is  restored  to  be  equal 
to  A  by  ATOB 
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10 


DO  10  J  -  1,  N 
A(J)  *  0. 

B(J)  -  0. 

DO  20  IND  «  1,  N0 
DO  20  J  *  1,  N 
II  -  J 
J1  -  J 

IF  (CABS  (A(J)).LT.  1.E-10)  THEN 
B(J)  -  CINC  (J) 

DO  40  I  «  1,M 
CALL  MAPINC 
40  CONTINUE 

CALL  LIKCHG 

.DLIK  »  -ALPHA  +  DLIK  *.5 
ELSE 

B(J)  -  0. 

DO  50  I  *  1,M 
CALL  MAPINC 
50  CONTINUE 

CALL  LIKCHG 

DLIK  *  ALPHA  +  DLIK  *.5 
END  IF 

IF  (DLIK  .GT.0)  THEN 
CALL  ATOB 
ELSE 

CALL  BTOA 
END  IF 
20  CONTINUE 

SUBROUTINE  MAPINC 

II  =  MAX{ 1 ,  II  -  L  +  1) 

J1  =  MIN(N ,  J1  +  L  -  1) 

DO  10  I  =  II,  J1 

IF  (CABS(B(I)),  GT.  1.E  -  10)  B(I)  =  CINC(I) 
10  CONTINUE 

RETURN 
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FUNCTION  CINC  (I) 

COMPLEX  CINC 

13  -  MAX  ( 1  ,  I  -  L  +  1 ) 

J3  -  MIN  (N,  I  +  L  -  1 ) 

B(I)  -  0. 

IP1  •>•1  +  1 
IM1  -  I  -  1 
DO  10  JI  *  13,  IM1 

10  B(I)  -  B(I)  -  C0NJG  <Q(I  -  JI  +  1))  *  B(JI) 

DO  20  JI  -  IP1,  J3 

20  B(I)  -  B ( I )  -  Q( JI  -  I  +  1)  *  B( JI ) 

CINC  «  (RMU<I)  +  B(X) )/Q( 1 ) 

RETURN 

SUBROUTINE  LIKCHG 

12  *  MAX( 1 ,  II  -  L  +  1) 

J2  »  MIN(N,  JI  +  L  -  1) 

DLIK  =  0 

DO  10  I  *  12,  J2 

10  DLIK  «  DLIK  +  REAL  (CONJG  (B(J)  -  A( J) )  *  RMU(J)) 

RETURN 

SUBFOUTINE  ATOB 
DO  10  I  -  1,  N 
10  B(I)  «  A(I ) 

RETURN 

SUBROUTINE  BTOA 
DO  10  I  =  1,  N 
10  A  (I )  =  B  (I ) 

RETURN 
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Appendix  B:  Derivation  of  the  steady  state  error  for  the  beet  linear 
estimate.  The  processes  {e^  :  j  e  z}  and  {a^  j  j  e  s)  (here,  X  denotes 
the  integers,  positive  and  negative)  have  spectral  representations  of  the  fora 

e.  «  J  eiXj  dZ(A), 

3  -IT 

a  «  /  eiXj  dA(  A) , 

3  -IT 

where  dz  and  dA  are  random  spectral  measures  satisfying  the  operational 
conventions 

E  dz( A)  dz(y)*  ■  ( a2/*)  6(A-y)dA, 

E  dA(  A)  dA(y)*  *  (pp2/*)  6(  A-y)  dA, 

E  dA(  A)  dz(  y)  -  0 

where  6(»)  denotes  Dirac's  delta  function.  See  section  3.3  and  equation 
(2.32)  of  Koopmans  (1974),  and  recall  that  {e^}  and  (a j }  are  complex 
valued  white  noises,  which  accounts  for  a  factor  of  2  in  the  latter 
equations.  The  BLE,  which  is  the  minimizer  of  (2.9),  is  given  by 

£  =  [bl  +  W*W]~1W*^ 

where 

(Bl)  b  =  o2/(pp2) . 

To  write  this  in  spectral  form,  we  first  define 

tf(  A)  =  £  w(k)e  iXk, 
k 


-34- 


where  w  is  the  waveform  of  (2.2).  Note  that  operating  with  W  is  just 
convolution  with  w.  For  the  particular  wave  for  used  in  the  simulations  of 
Section  5, 

L-1 

(B2)  0(A)  -  l  e~XAK 

k-0 


■  exp(-i(L  -  1)1/2]  sin  [LA/2]/sin[ A/2] . 


Now  it  is  easy  to  see  that 


[W4]  j  *  [  (bl  +W*W)_1  WWjr] 


J  — - -  eiAj  {d2  ( A)  +  v? (  A)  dA(  A)  }. 

-if  b  +  |«(A)r 


Here,  we  are  using  [W4] ^  to  denote  the  jth  component  of  W4.  The  steady 
state  error  (per  observation)  in  the  estimate  of  the  signal  Wa  is  then 


El  [W(4  -  a)]  |2  *  E|  /  -  I*—1  ■■  -  dZ(  A)  +  /  J±*i±L  dA(  A)  | 2  = 


-if  b  +  |0(A)  I ' 


-ir  b  +  10(  X)  | 


1  J  — gllHAlLi—  +  -£p..^l.a(  A)^2  dX__l  J  -mU.t  dA, 

-it  (b+|«(A)r)Z  (b+|*(X)|V  *  -ir  b  +  |v?(  A)  | 


i  /  - 


where  (Bl )  was  used  at  the  last  step.  Multiplication  of  this  latter  quantity 
by  N  gives  a  good  approximation  for  the  integrated  mean  squared  signal 
estimation  error  for  the  BLE  from  a  set  of  N  observations.  It  is  a  lower 
bound  since  the  boundary  effects  of  a  finite  observation  interval  always 
increase  the  MSE. 
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Appendix  Ct  Derivation  of  the  limiting  IMSE  for  the  MAP  estimate.  The 
integrated  mean  squared  error  (IMSE)  for  the  MAP  estimate  of  the  total  signal 
process  is  given  by 

N  2 

IMSE  -  E  l  | [W(A  -  a)]  |  , 

k*1 

where  A  denotes  the  MAP  estimate  of  a,  and  [Wu].  denotes  the  kc 

2  2 

component  of  the  N-vector  Wu.  Now,  we  will  fix  a  and  c  *  pp  ,  and  let 
p  +  0.  Under  these  conditions,  the  events  will  occur  less  frequently  and  the 
individual  signals  become  more  pronounced  (since  p  mast  tend  to  «•  ). 

Hence,  the  following  two  assumptions  will  become  approximately  valid: 

(i)  there  is  no  overlapping  of  signals; 

(li)  the  event  time  set  K  will  be  estimated  perfectly. 

In  any  case,  the  IMSE  derived  under  these  two  assumptions  will  be  a  lower 
bound  to  the  true  IMSE  for  any  p  >  0,  since  overlapping  signals  and 
uncertainty  in  K  will  only  increase  the  error. 

If  the  waveform  w  is  given  by  (5.1),  then  it  follows  from  (2.8)  and 
(2.7)  that  for  j  e  , 

3+lT1  ,  o 

4.  “  l  y./[(a"/p“)  +  l]  . 

3  k+j 

Hence,  the  IMSE  resulting  frc:.t  this  individual  signal  is 

L  E|A.  -  a. |2  =  L  [1  -  1/(o2/p2  +  L)]2  E|a  j  | 2  + 


=  2L02/(02/p2  +  L). 
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Since  the  expected  number  of  events  is  Np,  the  total  IMSE  ia  given  by 
IMSE  -  2NPLO2/(02/p2  +  L). 
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